%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function evaluates a layer potential with the combined field kernel.
% xx     is the matrix of target points, 
% C      is the contour.
% sigma  is the source distribution on C.
% kh     Helmholtz parameter.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function vv = evalpot(xx, C, sigma, kh)

n = size(C,2);
m = size(xx,2);

% Y_g1(i,j) = x_C(j) , X_g1(i,j) = x_eval(i)
% Y_g2(i,j) = y_C(j) , X_g2(i,j) = y_eval(i)
% Y_dg1(i,j) = dx_C(j)
% Y_dg2(i,j) = dy_C(j)

X_g1  = xx(1,:)'*ones(1,n);
X_g2  = xx(2,:)'*ones(1,n);
Y_g1  = ones(m,1)*C(1,:);
Y_g2  = ones(m,1)*C(4,:);
Y_dg1 = ones(m,1)*C(2,:);
Y_dg2 = ones(m,1)*C(5,:);

ima   = sqrt(-1);
% tangent vector : t = (dx,dy)
% -> normal vector n = (dy,-dx) / |t| 
% nn1(i,j) = nx(j)
nn1   = ( Y_dg2./sqrt(Y_dg1.*Y_dg1 + Y_dg2.*Y_dg2));
% nn2(i,j) = ny(j)
nn2   = (-Y_dg1./sqrt(Y_dg1.*Y_dg1 + Y_dg2.*Y_dg2));
% dd(i,j) = |x_C(j) - x_eval(i)|
dd    = sqrt((Y_g1 - X_g1).^2 + (Y_g2 - X_g2).^2);

% EVAL(i,j) = h |dx/dt(j)| (
%              -k*n(j).(x_C(j) - x_eval(i)) / d(i,j)*H1(k*d(i,j))
%              + i*k*h*H1(k*d(i,j)))

EVAL  = ((nn1.*(Y_g1 - X_g1) + nn2.*(Y_g2 - X_g2)).*(1./dd).*(-kh*besselh(1, kh*dd)) + ...
        ima*kh*besselh(0, kh*dd)).*sqrt(Y_dg1.^2 + Y_dg2.^2);
% vv(i) = sim_j EVAL(i,j)*sigma(j)
vv    = EVAL*sigma;

return
